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ABSTRACT 

Wc perform a scries of high-resolution N-body simulations of cosmological structure 
formation starting from Gaussian and non-Gaussian initial conditions. We adopt the 
best-fitting cosmological parameters from the third- and fifth-year data releases of 
the Wilkinson Microweive Anisotiopy Probe and we consider non-Gaussianity of the 
local type parameterised by eight different values of the non-linearity parameter /nl • 
Building upon previous work based on the Gaussian case, we show that, when ex- 
pressed in terms of suitable variables, the mass function of friends-of-fricnds haloes is 
approximately universal (i.e. independent of redshift, cosmology, and matter transfer 
function) to good precision (nearly 10 per cent) also in non-Gaussian scenarios. We 
provide fitting formulae for the high-mass end [M > 10^'^ h~^MQ) of the universal 
mass function in terms of /nl, and we also present a non-universal fit in terms of both 
/nl and z to be used for applications requiring higher accuracy. For Gaussian initial 
conditions, we extend our fit to a wider range of halo masses (M > 2.4 x 10^" h'^Mo) 
and we also provide a consistent fit of the linear halo bias. We show that, for realistic 
values of /nl, the matter power-spectrum in non-Gaussian cosmologies departs from 
the Gaussian one by up to two per cent on the scales where the baryonic-oscillation 
features are imprinted on the two-point statistics. Finally, using both the halo power 
spectrum and the halo-matter cross spectrum, we confirm the strong /c-dependence 
of the halo bias on large scales (fc < 0.05 /i Mpc~^) which was already detected in 
previous studies. However, we find that commonly used parametcrisations based on 
the peak-background split do not provide an accurate description of our simulations 
which present extra dependencies on the wavenumber, the non-linearity parameter 
and, possibly, the clustering strength. We provide an accurate fit of the simulation 
data that can be used as a benchmark for future determinations of /nl with galaxy 
surveys. 

Key words: cosmology: theory, dark matter, large-scale structure - methods: N-body 
simulations - galaxies: haloes, clusters. 



1 INTRODUCTION 

The detection of temperature anisotropies in the Cosmic Mi- 
crowave Background (CMB) provided evidence that large- 
scale structure formation in the universe was seeded by small 
density fluctuations generated at early times. The statisti- 
cal properties of these seeds are usually modelled with a 
Gaussian random field. Historically the Gaussian approxi- 
mation was introduced for mathematical convenience. In the 
absence of a solid model for the generation of density fluc- 
tuations the Gaussian hypothesis was accepted on the b asis 
of the central limit theorem (e.g. iBardeen et al.|[l986l and 



references therein). The advent of inflationary models pro- 
vided further support for Gaussianity. Small-amplitude cur- 
vature perturbations generated during a standard inflation- 
ary phase (singl e field, slow roll) ar e very nearly Gaussian 
distributed (e.g. iBartolo et al]|2004l and references therein). 

However, many variants of the inflationary scenario pre- 
dict appreciable levels of primordial non-Gaussianity. In 
terms of Bardeen's gauge-invariant potential, most o f 
these models (but not all, see e.g. ICreminelli et al.l |2007| ) 
can be reduced to the form: 

$ = '^ + /nl(</>'-(<^')), (1) 
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where (j> is an auxiliary Gaussian random field and /nl quan- 
tifies the amount of primordial non-Gaussianity. On sub- 
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horizon scales, "I> = —4' where ^ denotes the usual pecu- 
liar gravitational potential related to density fluctuations 
via Poisson's equation. The parameter /nl thus has the 
same sign as the skewness of the density probability distribu- 
tion function. This local form of non-Gaussianity (note that 
equation ((T)) applies in configuration space) can be obtained 
from a truncated expansion of the effective inflaton poten- 
tial (ISalopek fc Bondl Il990l : iFalk et al.l Il993l : iGangui et al] 
1 19941 ). The parameter /nl thus encodes information about 
the inflaton physics. Standard infiation gives |/nl| <C 1 
jSalopek fc Bondlll990l : iMaldacenalliooj ). However, even in 
this case, the non-linear evolution of perturbations on su- 
perhorizon scales yields an observable /nl of order unity 
(which, in reality, s hould be scale and redshift dep endent; 
iBartolo et all |2005| . see also iPvne fc Carrolll [l99g ). Large 
valu es of |/nl I natural ly arise i n multi-field inflation models 



(e.g. Linde fc Muk ha novlll997l : for an extensive review see 



iBartolo et aiTT2004l ') and even in cy clic or ekpyrotic mod- 
els of the universe with no inflation dCreminelli fc Senatord 
I2OO7I : iBuchbinder et al.ll2008l : iLehners fc Steinhardtll2008l ') . 

Observational constraints on /nl have been derived 
studying t hree-point statistics of tem perature fluctuations in 
the CMB (|Komatsu fc Spergelll200ll ). The recent 5-year data 
from the Wilkinson Microwave Anisotropy Probe (WMAP) 
jive —9 < ./nl < 11 1 at the 95 per cent confidence level 
Komatsu et al.|[2008l ). Parallel studies on the same dataset 
give —178 < /nl < 64 using Minkowski functionals (Ko- 
matsu et al. 2 008) and —8 < / nl < 111 from wavelet 
decomposition (|Curto et al.ll200a ). Some recent reanalyses 
of earlier 3-year WMAP data claim substantial evidence 
for positive /nl: 27 < /nl < 147 from the bispec trum 
of temperature fluctuations (|Yadav fc Wandeltl l2008h and 
23 < ./nl < 75 from th eir one-point distribution func- 
tion jjeong fc Smootll2007l ). On the other hand, a study of 
Mi nkowski functionals on the 3- year data gives — 70 < /nl < 
91 iHikage et al. I I2OO8I) . Higher quality data are needed to 
improve these constraints. The upcoming Planck satellite 
should be able to reduce the uncertainty in /nl to ~ 5 
(|Komatsu fc Spergel|[200ll ). 

Alternatively, one might use observational signatures 
of primordial non-Gaussianity imp rinted in the large-scal e 
structure (LSS) of the universe (e.g. lMoscardini et al.[l99ll ). 
Ideally, one would like to use high-redshift probes as the 
non-linear growth of density fluctuations quickly super- 
imposes a strong non-Gaussian signal onto the primor- 
dial one so that the latter might then be difficult to re- 
cover. For instance, the large-scale distribution of neutral 
hydrogen in the era between hydrogen reco mbination and 
reionisation encodes information on /nl (|Pillepich et al.l 
I2OO7I ). This could be probed by detecting the redshifted 
hyperfine 21-cm transition with very low-frequency ra- 
dio arrays from space. In principle, an experiment of 
this kind can limit /nl to A/nl < 1 (jPillepich et al.l 
I2OO7I . see also ICooravl 20061 ). However, it is not clear yet 
whether such an experiment will ever be possible due to 
technical complexity and problematic foreground subtrac- 
tion. At lower redshifts, /nl can be constrained prob- 
ing the statistics of rare events, as li ke as the mass 
funct i on of galaxy groups and clusters (iMatarrese et al] 
19861. I2OOOI: iKovama et all Il999l: iRobinson fc Bakeij bOOoF 



inconclusive due to low-numb er statistics (see e.g. IWillickl 
l2000l :[ Amara fc Refregieij|2004l and references therein). Even 
though cluster-mass estimates are still rather uncertain and 
massive objects are very rare, the observational perspectives 
look very promising. A number of galaxy surveys encom- 
passing large fractions of the observable universe are being 
planned (e.g. ground-based surveys as DES, PanSTARRS, 
and LSST, and the satellite missions EUCLID and ADEPT) 
and could potent i ally lead to solid mea surements of /nl (e.g. 
iDalal et al.ll2008l : ICarbone et"al]|2008l ). 

Primordial non-Gaussianity is also expected to modify 
the clustering properties of massive cosmic structures form- 
ing out of rare density fluctuations llGrinstein fc Wise! 19861: 
Matarrese et al. I Il986l : iLucchin et all Il988l : IKov' ama et al.l 
19991 ). Also in this case, however, the non-linear evolution 



Robinson et al.ll2000l : iLoVerde et~a]| I2OO8I ) . Early attempts 



of using cluster counts to constrain /nl have been rather 



of the mass density generally superimposes a stronger sig- 
nal than that generated by primordial non-Gaussianity onto 
the galaxy three-point statistics. The gal axy bispectrum is 
thus sensitive to /nl only at high redshift jVerde et al.|[2000l : 
IScoccimarro etaL II2OO4I: ISefusatti fc Komatsull2007l ). 

Recently, I Dalai et al.l (|2008l ) have shown analytically 
that primordial non-Gaussianity of the local type is ex- 
pected to generate a scale-dependent large-scale bias in the 
clustering properties of massive dark-matter haloes. This is 
a consequence of the fact that large and small-scale den- 
sity fluctuations are not independent w hen /nl 0. Simi- 
lar ca lc ulations have been p resented by Matarrese fc Verd3 
(l2008l).[siosar et al.l (|2008D . lAfshordi fc ToUevI (l2008l). and 
iMcDonaidI ( 20081 ) . Numerical simulations by iDalal et al.l 
( 20081 ) are in qualitative agreement with the analytical pre- 
dictions conflrming the presence of a scale-dependent bias. 
Using these analytical models for halo biasing to describe the 
clustering amplitude of luminous r ed galaxies and qua sars 
from the Sloan Digital Sky Survey, ISlosar et all (|2008l ) ob- 
tained —29 < /nl < 69 at the 95 per cent confidence level. 
This shows that LSS studies are competitive with CMB ex- 
periments to constrain primordial non-Gaussianity but also 
calls for more accurate parameterisations of the mass func- 
tion and clustering statistics of dark-matter haloes arising 
from non-Gaussian initial conditions. 

Most of the ana lytic derivations of the non-Gaussian 
halo mass function ( Matarrese et al ] I2OOOI : iLoVerde et ail 
|2008|. e.g.) are based on the extended Press-S chechter model 
( Press fc Schechteij [l973 : iBond et all Il99lh which, in the 
Gaussian case, is k nown to produce inaccurate estimates 
of ha lo abundance (|Sheth fc Tormeiil 1 19991 : Ijenkins et al.l 
I2OOII ). Similarly, the scale dependent bias is obtained ei- 
ther using the peak-background split model (jSlosar et al.l 
|2008| ) or assuming that haloes form from th e highest lin- 
ear density peaks (jMatarrese fc Verdel I2OO8I ). Both tech- 
niques have limited valid i ty in the Gaussian ca se (jjinel 
1 19981 : iPorciani et al. I ll999l : ISheth fc Tormenlll999l ). In this 
paper we test the accuracy of the excursion-set model 
and the peak-background split in the non-Gaussian case . 
This extends the previous studies of | Kang et al. I (|2007h . 
iGrossi et all (^2007^ andlDalal et al.l (|2008l ) for the halo mass 
function and of iOalal et al.l (|2008F for the halo bias by ex- 
ploring more realistic values for /nl with simulations of bet- 
ter quality. In practice, we run a series of high-resolution N- 
body simulations where we follow the process of structure 
formation starting from Gaussian and non-Gaussian initial 
conditions. The halo mass function and bias extracted from 
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Table 1. Specifics of tlie N-body simulations. 



Name 


/nl 


Npart 


Lbox 




Lsoft 


^start 


Cosmology 








(h-iMpc) 




(/i-^kpc) 






i. / OU 


+ ( OU 




izUU 


1 o/ifi 1 nil 
1.z4d X iU 


on 


OU 


WiVlAr^O 


1.500 


-f500 


1024^ 


1200 


1.246 X 10" 


20 


50 


WMAP5 


1.250 


+250 


1024^ 


1200 


1.246 X lO^i 


20 


50 


WMAP5 


1.80 


+80 


1024^ 


1200 


1.246 X lO^i 


20 


50 


WMAP5 


1.27 


+27 


10243 


1200 


1.246 X lOii 


20 


50 


WMAP5 


1.0 





10243 


1200 


1.246 X 10" 


20 


50 


WMAP5 


1.-27 


-27 


10243 


1200 


1.246 X 10" 


20 


50 


WMAP5 


1.-80 


-80 


10243 


1200 


1.246 X lOii 


20 


50 


WMAP5 


2.0 





10243 


1200 


1.072 X 10^1 


20 


50 


WMAP3 


2.750 


+750 


10243 


1200 


1.072 X IQii 


20 


50 


WMAP3 


3.0 





10243 


150 


2.433 X 10* 


3 


70 


WMAP5 


3.250 


+250 


10243 


150 


2.433 X 10** 


3 


70 


WMAP5 



Table 2. Assumed cosmological parameters. 



Name 


h 












WMAP3 
WMAP5 


0.73 
0.701 


0.76 
0.817 


0.95 
0.96 


0.24 
0.279 


0.042 
0.0462 


0.76 
0.721 



the simulations are then compared with the existing ana- 
lytical models and used to build accurate fitting formulae. 
These will provide a benchmark for future determinations 
of non-Gaussianity with galaxy surveys. 

The paper is organised as follows. In Section 2 we de- 
scribe our N-body simulations. In Sections 3, 4, and 5 we 
present our results for the halo mass function, the matter 
power spectrum and the halo bias, respectively. In Section 
6 we discuss th e impl ications of our results for the analysis 
bv lSlosar et al.l (|2008l ). Our conclusions are summarised in 
Section 7. 



2 N-BODY SIMULATIONS 
2.1 Specifics of the simulations 

We use the le an version of the tree-PM code Gadget-2 
l|Springelll2005l l kindly made available by Volker Springel to 
follow the formation of cosmic structure in a flat ACDM 
cosmology. We run three diflerent series of simulations 
(each containing 1024^ coUisionless particles) that differ in 
the adopted cosmology, box size (and thus force softening 
length, Lsoft), and initial redshift (details are summarised in 
Table [T| . The assumed cosmological parameters are listed 
in Table [2] For our series #1 and #3 they coincide with 



the 5-yr WMAP best estimates dKomatsu et al.ll2008l). The 
combined 3-yr WMAP+LSS results bv lSoergel et air(|2007l l 
are instead used for series #2. 

We produce non-Gaussian initial conditions di- 
rectly applying equation ([!]) after having generated the 
Gaussian random field <p with standard Fourier tech- 
niques. We consider eight values for the parameter 
/nl: -80, -27, 0, +27, +80, +250, +500, +750. The first five 
are within the cur rent constraints from CMB data 
(|Komatsu et al.l I2OO8I ) , while the three largest values are 



useful to compare with previous work. Within each series of 
simulations, we use the same set of random phases to gener- 
ate the Gaussian potential 0. This facilitates the comparison 
between different runs by minimising sample variance. 

The linear matter tr ansfer function, T( k), is computed 
using the Linger code (|Bertschingeij [2OO1I ) and is applied 
after creating the non-Gaussian potential $ in equation 
Particle displacements and velocities at Zstart are gener - 
ated using the Zel'dovich approximation (|Zerdovichlll97Gl ). 
A critical discussion of this choice is presented in the Ap- 
pendix. 

Particle positions and velocities are saved for 30 time 
steps logarithmically spaced in {1 + z)~^ between z = 10 
and z = 0. Dark-matter haloes are identified using the 
standard friends-of-friends (FOF) algorithm with a linking 
length equal to 0.2 times the mean interparticle distance. 
We only considered haloes containing at least 100 particles. 

Our first two series of simulations only include large 
periodic boxes covering a volume of (1200 /i~^Mpc)'' where 
we can study haloes with masses ranging from 10^'' up to 
10^^ h~^MQ. These simulations will be used to analyse both 
the mass function and the bias of dark-matter haloes. On 
the other hand, the third series includes simulations covering 
a volume of (150 h~^Mpc)^ . They will be used to study the 
mass function and the bias of low-mass haloes with lO'^" < 
M < 10" h'^MQ. 



2.2 A note on the definition of /nl 

The definition of /nl given in equation ^ depends on the 
cosmic epoch at which it is applied. The reason for this 
time dependence is that both potentials $ and (f) decay with 
time proportionally to g{a) = D{a)/a with D{a) the linear 
growth factor of density fiuctuations and a the Robertson- 
Walker scale factor. 

In this paper, we define /nl by applying equation ([T]) at 
early times, namely at z = 00. Oth e r aut hors have adopted 
different conventions. iGrossi et al. I (|2007h use the linearly- 
extrapolated fields at z = to define /nl. Therefore, their 
values of the /nl parameter need to be divided by the fac- 
tor g{oo) / g{0) to match ours. In the WMAP 5 cosmology , 
g{oo)/g{0) ~ 1.3064. On the other hand. lDalal et al.l (|2008l ) 
apply equation ([1} at Zstart, the redshift at which they gen- 
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Figure 1. The universal mass function in our Gaussian simula- 
tions Run 1.0 (triangles), Run 2.0 (squares), and Run 3.0 (circles) 
is compared with a number of fitting formulae listed in Table [S] 
Data are equispaced in Incr"^ and only bins containing more than 
30 haloes are shown. The vertical do t ted li n es indicate t he up - 
per mass limits used in Jenkins et al ] l|200lh . lReed et al.l l|2003l ). 



and lWarren et"aL I i l2oodr The corresponding low-mass limits are 
all equal or smaller than Incr"! = —1.2. The lower panel shows 

between our data points and the 



the residuals 



(data— fit) 
data 



different fitting functions. Here we only show data with a Pois- 
son uncertainty better than 5 per cent. For clarity only outputs 
from Run 1.0 (triangles, Incr^^ > —0.3) and Run 3.0 (circles, 
ln(T~l < —0.3 ) are plotted. 



erate the initial conditions for the simulations. This agrees 
with our definition to better than 0.01 per cent. 

The sign convention for the non-linearity parameter 
might possibly generate further ambiguity. In our simula- 
tions, positive values /nl correspond to positive skewness 
of the mass-density probability distri bution function. Th e 
same convention has b een adopted by Grossi et al.l |20o3), 
iKang et all (|2007l ) and lDalal et al] (|2008l ). 



3 THE HALO MASS FUNCTION 

One of the long standing efforts in cosmology is to determine 
the mass function of dark matter haloes dn/dM{M, z) ~ i.e. 
the number of haloes per unit volume per unit mass at red- 
shift z - from the statistical properties of the linear density 
field. Analytical work has suggested that, when expressed 
in terms of suitable variables, the functional form of dn/dM 
should be universal to changes in redshift and cosmology 
JPress fc Schechteill97i : lBond et al.ll99ll : [Sheth fc TormenI 
1 19991 ). N-body simulations have shown that this is approx- 



imately true when structure formation is seeded by Gaus- 
sian perturbations (Ijenkins et aL 2001 : Evrard et ahlbood : 
IWhite|[200^ : IWarren et al.ll2006l : iTinker et al.ll2008l ). 

Following these studies, we describe the halo abundance 
in our simulations through the following functional form 



d \n[a-^{M,z)] 
dM 



(2) 



where pm is the mean background matter density today, and 
(T^(M, z) is the variance of the linear density field 

1 r^^ 

a\M, z)^— Pnn{k, z) W^{k, M) dk, (3) 

^1" Jo 

with Pun{k,z) the corresponding power spectrum and 
W^{k,M) some window function with mass resolution M 
(here top-hat in real space). The validity of equation ([2]) 
has been widely tested against numerical simulations and 
useful parameterisations for /(a) ha ve been provided start- 
ing from Gaussian initial conditions ((Slieth fc Tormen|[T999l : 
Ijenkins et al]|200ll : IWarren et al.ll2006l ). These fitting func- 
tions have an accuracy ranging from 5 to 20 per cent de- 
pending on redshift, c osmology, and the exact definition of 
halo masses. Recently. [Tinker et al.l (|2008l ) have detected de- 
viations from universality in /(cr): redshift-dependent cor- 
rections are needed to match the mass function in simula- 
tions with an accuracy of 5 per cent. This result is based on 
haloes identified with the spherical overdensity algorithm. It 
is well known that the mass function of EOF haloes shows a 
more universal scaling even though other halo finders might 
be more directly linked t o actual observables (| Jenkins et al.l 
I2OOII : [ Tinker et al.ll20o3 ). Deviations from universality for 
FOF haloes will be further discussed in Section 3.3. One 
should anyway keep in mind that baryonic physics can cause 
30 per cent devia tions in dn/dM wi th respect to the pure 
dark-matter case (jStanek et al.|[2008l ). 



3.1 Halo mass function from Gaussian initial 
conditions 

The halo mass functions extracted from our Gaussian sim- 
ulations - Run 1.0 (triangles). Run 2.0 (squares), and Run 
3.0 (circles) - are presented in Figure [T] The combination 
of difi'erent box sizes allows us to cover the very wide range 
— 1.2 < lna~^ < 1.1 which roughly corresponds to the mass 
interval 2 x 10^° < Af < 5 x 10^^ h'^M^ sX z = Q. Figure[T] 
has been obtained by combining data from snapshots at red- 
shifts 2 < 1.6. Note that, at a fixed redshift, larger values of 
correspond to higher masses. On the other hand, with 
increasing the redshift, larger values oi are associated 
with a given halo mass. Even though datapoints correspond 
to difi'erent redshifts and cosmologies, they all form a well 
defined sequence. This indicates that the function /(cr) is 
universal to good approximation. For a given cr, outputs at a 
fixed redshift scatter around the universal sequence by 10-15 
per cent. A number of fitting formulae have been proposed 
in the literature to parameterise this sequence. In Figure [TJ 
we compare some of them (summarised in Table [3| with our 
datapoints. Fractional deviations between models and data 
are shown in the bottom panel. Barring the classical Press- 
Schechter result, all the fitting formulae describe our data 
to better than 20 per cent . The best agreemen t is found all 
over the mass range with IWarren et al.l (|2006l ) foUowed by 
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Table 3. Widely used parameterisations for the halo mass function deriving from Gaussian initial conditions. 
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Acronym 


Reference 


PS 


Press & Schechter (1974) 


ST 


Sheth & Tormen aggQ) 


J 


Jenkins et al. (2001) 


R 


Reed et al. (2003) 


W 


Warren et al. (2006) 


T 


Tinker et al. (2008) 



Functional form Parameters 



/PSM = ^ exp i^-^j = 1.686 

/sTM = Ayf ^^exp(^-^) [l+(^^Y^ A = 0.322, a = 0.707, p = 0.3 

/,(o-) = Acxp (-|lno-i +B|P) A = 0.315, B = 0.61, p = 3.8 

/RH = /sT{-)exp(^^j^^^) a = 0.7,6 = 5 

/w(o-) = A {a-'' + b) cxp (-^) A=0.7234, a=1.625, b=0.2538, c=1.1982 
f'Y(cr) = A l^(^) " + ij exp ^— vary with halo overdensity 



Ijenkins et al.l (|200ll ) which both s how deviations from our 
data at the 10 per cent level. The ISheth fc TormenI (|l999l ') 
model also provides an accurate description of the data for 
small halo masses but tends to overestimate the abundance 
of the mo st massive objects. On the other hand, the fit by 
iReed et~al . (2003.) tends to underestimate the high-mass tail 
of the mas s function. Overa l l our findi ngs are in good agre e- 
ment with iHeitmann et al. I (|2006D and lLukic etahl (120071) 

FoUowing lWarren et al.] (|2006l ) and lTinker et all (|2008l) . 
we fit the outcome of the simulations with the function 



tt-tf-^ 



D + B ( - 

a 



(4) 



The best-fitting parameters have been determined through 
minimisation using the Markov Chain Monte Carlo 
method, and read: 

A = 1.868 ±0.019 

B = 0.6853 ± 0.0035 (5) 
C = 1.2266 ± 0.0049 
D = 0.2279 ± 0.0022 . 

In ter ms of the parameterisation given in IWarren et al.l 
l|2006l ) and reported in Table O this corresponds to 
(A, a, b, c) = (0.6853, 1.868, 0.3324, 1.2266). The fit in equa- 
tion ^ describes our dataset up to deviations of a few per 
cent over the entire mass and redshift ranges for Run 1.0 
and Run 2.0, while it shows larger deviations (up to nearly 
10 per cent) towards the high-mass end of Run 3.0, (see Fig- 
ure [TJ. It is important to remember, however, that Run 3.0 
covers a much smaller volume than the others and thus is 
more severely affected by sample variance. 



3.2 The universal halo mass function from 
non-Gaussian initial conditions 

Is the function f{a~^) universal also in the non-Gaussian 
case? This question is addressed in Figure [2] where we show 
the output of our main series of simulations at four redshifts 
{z = 0,0.5, 1, 1.6) to test the scaling of the mass function 
in terms of a~^. Only bins containing at least 20 haloes 
are considered. Within a certain tolerance, the halo mass 
functions at different masses and redshifts all lie on the same 
curve for a given /nl- The scatter of the points at a fixed 
redshift around this curve roughly amounts to 10 per cent, 
and it becomes smaller towards our largest values of /nl. 
We thus generalise equation ((2)1 to non-Gaussian initial 
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Figure 2. Universality of the mass function arising from non- 
Gaussian initial conditions. Colors refer to simulations with dif- 
ferent values of /nl as indicated by the labels. Symbols identify 
the redshift of the simulation output from which the mass func- 
tion has been calculated, namely z = (triangles), 0.5 (circles), 
1 (squares), 1.6 (diamonds). 



Pm 



d ln[(T-^(M,z)] 



(6) 



conditions by assuming that 

^(/..,M,.) = /(/N.,.)^^ 

and we provide a fitting formula for /(/nl,'^). Given the 
similarity to the Gaussian case, we still adopt the functional 
form given in equation but let the parameters A, B, C, D 
vary with /nl- The best-fitting values have been determined 
in two steps. We first used a Markov Chain Monte Carlo 
method to determine A, B, C, D at fixed /nl through 
minimisation. The results suggest that the /nl dependence 
for each parameter of the mass function can be accurately 
described by polynomials of different orders. Eventually, we 
used the data to derive the coefficients of these polynomials. 



The degree of complexity required to fit the simulation 
data grows considerably with increasing /nl. For —80 ^ 
/nl 250 (a range that fully encloses the values currently 
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Figure 3. Comparison between the halo mass function from our main series of simulations (triangles) and the corresponding fitting 
functions (lines). Values —80 /nl ^ 250 and the fit in equation l[7} are considered in the left panel. All the simulations and the 



polynomial fit in equations (|8} and l(9} are shown in the right panel. The lower panels show residuals 
with a statistical uncertainty which is smaller than 5 per cent. 



(data — fit) 
data 



for data points 



Parameter pi p2 



A 


1.694 


-0.00199 


B 


0.566 


-0.00029 


C 


1.151 


-0.00071 


D 


0.287 


-0.00030 



Table 4. Best-fitting values for the linear coefficients of the 
universal mass- function parameters given in equation 0. The 
quoted values are truncated at the first digit which is affected by 
the statistical errors. This provides an accurate description of our 
simulations for —80 /nl ^ 250. 



Parameter 


PI 


P2 


P3 


P4 


P5 






(10-3) 


(10-7) 


(10-9) 


(10-12) 


A 


1.708 


-2.07 


3.1 








B 


0.560 


-0.46 


-f 12.46 


-2.36 


-f2.65 


C 


1.150 


-0.76 


-1-2.7 








D 


0.293 


-0.16 


-14.07 


+3.88 


-3.94 



Table 5. As in Table |4] but for the fitting formula in equations 
||8} and This accurately describes the mass function in all our 
non-Gaussian simulations (—80 ^ /nl ^ 750). 



allowed by CMB studies), the mass-function parameters in 
equation Q are well approximated by the linear relation 

P(fNL) =pi +P2 -fNL, for P = A,B,C,D. (7) 

Table U lists the corresponding best-fitting parameters. 
The quality of this fitting formula is assessed in the left 
panel of Figure [S] where the mass function for the simula- 
tions with /nl = -80, -27, 0, +27, +80, +250 is compared 
with the corresponding fit. Residuals (shown in the bot- 
tom panel) are smaller than 5 per cent all over the range 
—0.2 < Inff-i < 0.8 corresponding to the mass interval 
2 X 10^3 < A/ < 2 X IQ^^h'^MQ at 2 = 0. 
On the other hand, equation ((7} is not suitable to account 
for values of /nl substantially larger than 250. To obtain 
an accurate fit of the universal halo mass function over the 
range — 80 ^ /nl 5; 750 we had to consider polynomials up 



to 4**" order in /nl: 

P(fNL)=Pl+P2-fNL+P3-fNL, for P = A, C (8) 

and 

P(fNL) = Pi + P2 ■ fNL + P3 ■ fNL + P4 ■ fNL + PS ■ fNL, (9) 

for P = B,D. 

The best-fitting values of the parameters above are listed 
in Table [5] while the corresponding functions are compared 
with the simulation data in the right panel of Figure [S] 
Also in this case residuals are smaller than 5 per cent for 
In cr-^ < 0.8. 

The universality of the fitting formula in equation © 
has been further tested against our non-Gaussian simula- 
tion of the WMAP3 cosmology, Run2.750, which has not 
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been used to determine the best-fitting parameters. This 
bhnd check shows that, in the range -0.27 < Ina-i < 
0.94 (roughly corresponding to 1.6 x 10^^ < M < 2.2 x 
10^^ h~^MQ at z = 0), the provided fit reproduces the mass 
function with an accuracy of 5 per cent. 

We warn the readers against extending our fitting for- 
mulae beyond their range of validity, in particular at low 
halo masses. The simulations of our main series resolve 
10^^ ft-^Me haloes with 100 particles. For /nl / 0, our 
analytical formulae for the mass function have been derived 
using only haloes that are more massive than this limit. 
Moreover, since the high-mass tail of the mass function is 
enhanced (suppressed) for positive (negative) values of /nl 
with respect to the Gaussian case, mass conservation re- 
quires that the opposite effect is seen at lower masses. We 
have directly tested the goodness of our fit towards the 
smaller masses using Run3.250 (which has a boxsize 8 times 
smaller than for the simulations in the main series but the 
same number of particles) and indeed found that the fitting 
formulae in equations ((Sjl and ((9)1 systematically overesti- 
mate the abundance of small mass haloes by 10-30 per cent. 
We will address the low-mass tail of the mass function for 
/nl / in future work. 

On the other hand, for Gaussian initial conditions, we 
combined simulations with different box sizes to derive the 
fitting function in equations Q and ([5]). This allowed us to 
extend the validity of our fit to the much wider mass range 
2.4 X 10"' < A/ < 5 X 10^^ /i-^Mq. 

Our fitting formulae give three different approximations for 
the universal mass function in the Gaussian case. In general, 
the fit given in equations ^ and ((5)1 has to be preferred 
as it has been obtained from a richer dataset spanning a 
much wider range of halo masses. However for masses above 
lO^^/i^^M© at z = 0, the fit in equations ^ and (O and 
Table |4] provides the most accurate representation of our 
data. In any case, the different fitting functions never deviate 
by more than 3-4 per cent. Also note that our two fitting 
functions for the non-Gaussian simulations agree by better 
than 1 per cent for — 27 ^ /nl ^ 80 and by a few per cent 
for /nl = -80 and /nl = +250. 



3.3 The limit of universality: redshift dependence 

Regardless of the value of /nl, we have found that the halo 
mass function is universal, when written in terms of , 
with an accuracy of roughly 10 per cent. If one is interested 
in giving analytical approximations for the halo mass func- 
tion which are more accurate than the universal fit, it is 
nece ssary to introduce redshift-dependent corrections (see 
also [Tinker et al.l |2008| for the Gaussian case) . In the left 
panel of Figure |4l we show how well the universal fit (whose 
parameters are listed in Table [5]) describes the simulation 
outputs at z = 0, 0.5, 1, 1.61. At 2: = and for masses 
M ^ 4 — 5 • IO^^Mq the fitting formula deviates for the data 
by more than 10 per cent. The smaller the redshift, the worse 
is the agreement between the data points and the universal 
fit. The bigger the /nl, the less critical is the comparison. 
In this Section we provide a non-universal fit which is very 



Parameter 


Pi 


P2 


P3 


P4 






(10-1) 


(10-1) 


(10-4) 


A 


1.82 


2.85 


4.53 


-5.92 


B 


0.578 


5.30 


7.53 


-7.73 


C 


1.15 


9.52 


9.08 


-4.42 


D 


0.294 


4.92 


4.67 


-3.80 



Table 6. Best fitting values for the mass-function parameters 
given in equation l|ll|l . This provides an accurate description of 
the data for ^ z ^ 0.5 and —80 ^ /nl ^ 80. The universal 
function of equations l(6]l, l|4}, Jsj. and Q should be used other- 
wise. 



accurate at low redshift. In particular, we write: 

'^'^ ff M .\ f(f ^ , Pmdhi[oo^^(M)| 
— (/NL,M,.) = /(/NL,ao,z) — — , (10) 

where ctq = a{z = 0) = g[z) / Dj^{z) is the rms deviation 
of the linear density field at z = 0. We approximate / with 
the functional form given in equation (|4]) but now let the pa- 
rameters A, B, C, D vary with both /nl and z. Markov Chain 
Monte Carlo fitting suggests that each parameter A, B, C, D 
of the mass function can be accurately described as follows; 

P(z, fNL) = Pi [1 + P2 ■ Z -f P3 • Z^] [1 -I- P4 • Inl] • (H) 

The best fitting parameters for —80 ^ /nl ^ 80 and 
^ 2 ^ 0.5 are listed in Table |S] while the quality of 
the fitting formula is assessed in the right panel of Figure 
|4l Residuals are smaller than 5 per cent all over the mass 
range, indicating that for —80 ^ /nl ^ 80 and ^ z ^ 0.5 
the fit of equations 1101 U and [TT] has to be preferred to 
the universal fit given in the previous Section. On the other 
hand, for higher values of |/nl| and for higher redshifts, the 
universal fit gives a better and more economic (in terms of 
parameters) description of the data. 



3.4 Comparison with theoretical models 

The halo mass function arising from mildly non-Gaussian 
initial conditions can be modelled by generalizing the Press- 
Schechter formalism. Using the saddle-point approximation 
to evaluate the probability for the linear density field to be 
above a given threshold value, iMatarrese et al.| (I2OO0I) have 
derived a model for dn/dM. More recentlv. iLoVerde et all 
(|2008l ) presented another expression for the mass function 
by using the Edgeworth asymptotic expansion for the 
probability density function of the linear density field. In 
both cases, only leading-order corrections in /nl have been 
accounted for. In absolute terms, these models are not 
expected to be accurate as they should suffer from the same 
shortcomings as the Press-Schechter model in the Gaussian 
case. However, they can be used to compute the fractional 
non- Gaussi an correctio n /(M, z, fNh )/ .f(M, z, /nl = 0) 
(jVerde et al . 2000 : LoVerde et a]||2008l l. In Figured we use 
this quantity to test the models against our simulations. 
Datapoints with errorbars show the N-body output at 
z = 0,0.5 and 1, wh il e the dotted lines indicate t he mo dels 
of IMatarrese et all (|2000l ) and iLoVerde et"al (|2008l ') as 
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Figure 4. Left panel: Mass function residuals of Runl.-80, Runl.80, Runl.250, Runl.750 with respect to the universal fit given in 
equations (jsjl and ^ at redshift 2 = 0, 0.5, 1, 1.61 (indicated by the symbols and colors). Only data points with a statistical error 
smaller than 10 per cent are shown. Right panel: As in the left panel but for the non-universal fit given in equation (lilt . In this case, 
only data points with an accuracy better than 5 per cent are shown. 



indica ted in equations (B.6|3 and (4.19) of iLoVerde et all 
l|2008h . respectively. It is evident that the models overesti- 
m ate the non- Gaussian c orrection. Following the indications 
in ICarbone et al. I (|2008D , we also show a modified version 
of the models which is obtained by lowering the critical 
threshold for halo collapse as (5c — 1.5 (solid lines in Figure 
O. Such a correction vastly improves the agreement with 
the simulat i ons. 

iDalal et al.l (|2008l ) proposed to fit the halo mass function 
in terms of the convolution between dn/ dM{f^i^ = 0, M, z) 
and a Gaussian kernel in M with a /NL-dependent mean 
and variance. Figure [5] shows that their fit tends to over- 
estimate the non-Gaussian corrections especially for large, 
positive values of /nl and masses M < IQ^'^h^^ Mq. On 
the other hand, for |/nl| < 100 it has a similar accuracy as 
the formulae derived from the Press- Schechter formalism 
corrected with the reduced threshold. 

The good agreement between the fractional non- 
Gaussian corrections derived from the modified PS mod- 
els and from the simulations is not enough to derive /nl 
from future observations of galaxy clusters. In fact the ra- 
tio /(^, /nl)//(-z, /nl = 0) is not an observable: the only 
quantity that we can hope to compare with observations is 
the mass function. In order to make predictions for dn/dM, 
the models for the fractional non-Gaussian correction need 
to be multiplied by a Gaussian mass function. This step 



might introduce relatively large systematic errors (see Figure 
1) which could degrade any measurement of /nl based on 
the cluster mass function. We address this issue in Figure [6] 
where we plot the fractional deviation of some model predic- 
tions for the function / with respect to the simulation output 
(results are very si milar for different valu es of /nl)- We con- 
sider the model bv lLoVerde et aj (|2008l ) corrected with the 
factor J\f and mul t iplied by three differe n t Gau ssian models: 
ISheth fc TormenI (|l999l 'l. IWarren et all (|2006l '). and our fit 
with /nl = 0. Note that some of the final outcomes system- 
atically differ by 10-20 per cent over the entire mass range 
covered by the simulations. This clearly shows that a care- 
ful measurement of the Gaussian mass function is necessary 
to avoid a biased estimation of the non-line arity parameter. 
Note that, for |/nlI < 100, the models by iMatarrese et al.l 
(|2000l ') and ILo Verde et all (|2008l ) (both with the reduced col- 
lapse threshold) combined with our Gaussian fit are in rather 
good agreement with the numerical mass fun ctions (similar 
result s are obtained using the Gaussian fit by I Warren et al.l 
(|2006l ) for masses below a few xlO^''ft~^M0). Perhaps not 
surprisingly, no model describes the simulation data for all 
the values of /nl as well as our fitting formulae for the non- 
Gaussian mass function given in Sections 3.2 and 3.3. 



^ This fixes a typo in equation (68) of lMatarrese et al.l | |200CJ) 
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Figure 5. Compari son between tlie lialo mass functi ons from our simulations and from the models by iMatarrese et alj l|2000t l. by 
I Lo Verde et all l l2008h . and the fit by iDalal et alJ bOOSi ) for different values of /nl (different panels) and for z = 0, 0.5, 1 (triangles, 
ci rcles, squares, respecti yely) . The q uantity which is plo tted is the ratio /(z, /nl)//(2, /nl = 0, ). The dotted lin es indicate the model s 
of lMatarrese et al ]|2000l (green) and lLo Verde et alllioosl (magenta), as they appear in equations (B.6) and (4.19) of lLo Verde et a] | |2008|) , 
respectiyely. The corresponding solid lines indicate the same model s with a reduced th reshold for halo collapse: 5c ~ 1.5. The blue solid 
lin es are obtained by c onvolving the /NL-dcpcndcnt kernel given in lDalal et al. I ||2008|) with the mass-function fit for the Gaussian case 
bv lWarren et al.l l|2006h . 



3.5 Summary of accuracy and range of validity of 
the mass function fits 

In order to facilitate the use of our fitting formulae for the 
halo mass function we summarize here their accuracy and 
range of validity. 

• For -80 ^ /nl s; 80 and ^ z ^ 0.5 the best de- 
scription (with 5 per cent accuracy) of our numerical data 
is given by equations (fTO)) . ^ and (fTTj) : 

• For larger values of /nl and z (but with /nl ^ 750 and 
z ^ 1.6) or whenever an accuracy of 10 per cent is enough, 
the universal fits of Section 3.2 should be used: 



- universal fit for —80 ^ /nl ^ 250: equations (|4]), ((Tjl 
and Table H 

- universal fit for —80 < /nl ^ 750: equations (HI, (|8l), 
© and Tabled 



4 MATTER POWER SPECTRUM 

In this section we study how non-Gaussian initial conditions 
influence the power spectrum of the mass density field. At 
tree level, the power spectrum does not depend on /nl in 
Eulerian perturbation theory. However, one-loop corrections 
make the power spectrum /NL-dependent. Qualitatively, 
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Figure 6. Residuals between the simulated mass functions and various model prescriptions, for Runl.80 (left-hand panels) and Runl.500 
(right-hand panel) at 2 = 0, 0.5, 1. The models have been obtained multiplying the formula for /(z, fj^ i^)/ f{z, = 0; ) bylLo Verde et all 
1120081) (with the r educed Sc as the solid lines in Figure [Sj with different Gaussian mass functions: ISheth fc TormenlllQQgT (magenta) , 
IWarren et al ]|2006l (red), and our Gaussian fit (green). The black lines show residuals with respect to our fitting formula given in Sections 
3.2 and 3.3. Only data with errors smaller than 10 per cent arc shown. 



theoretical expectations are that positive (negative) values 
of /nl tend to enhance (suppress) the ampUtude of the 
power spectrum on non-Unear scales. In Figure[7]we plot the 
ratio of power spectra P(fc, /nl)/^'!^, /nl = 0) extracted 
from the simulations of our main series at redshifts z = 
and 1. The matter power spectrum of non-Gaussian models 
appears to deviate already by a few per cent at k — 0.1 h 
Mpc~^. As expected, deviations become more severe with 
increasing the wavenumber k. Our res ults are in agreemen t 
with the perturbative c alculations by pTaruva et al.l (|2008h . 
We note, however, that iGrossi et al. I (|2008t ) found smaller 
deviations between the non-Gaussian and Gaussian power 
spectra at larger values of k and /nl- 

Our results have two important practical implications. 
First, the widespread habit of using the Gaussian matter 
power spectrum to determine non-Gaussian bias parameters 
leads to scale-dependent systematic errors that might be- 
come severe when high-precision is required. Second, primor- 
dial non-Gaussianity modifies the power-spectrum on the 
scales where baryonic oscillations (BAOs) are present. Re- 
versing the argument, two-point statistics could be also used 
to constrain the value of /nl - Note however, that all probes 
based on galaxy clustering will suffer from uncertainties in 
the bicis parameter (and its scale dependence) that might 
hinder a measure of /nl based on the study of BAOs. On 
the other hand, weak lensing studies will directly measure 
the matter power spectrum. The target of many future wide- 
field missions is to provide estimates at the per cent level. 
For parameter estimation, a comparable accuracy will be 
required on model spectra within a wide range of wavenum- 



bers centred around fc ~ 1 ft Mpc ^ (|Huterer fc Takadal 
l2005l ). Therefore, even values of /nl within the current CMB 
constraints could imprint detectable effects in the matter 
power spectrum at the scales of interest. The key question is 
whether one can discern the effect of /nl and, consequently, 
how much primordial non-Gaussianity will affect the esti- 
mate of the other cosmological parameters. We will get back 
to this in future work. 



5 HALO CLUSTERING 

The clustering of dark-matter haloes is biased relative to 
that of the underlying mass distribution by an amount which 
depends on halo mass, redshift, and the scale at which 
the clustering is considered (see e.g . IMo fc White! Il996l : 
ICatelan et"aiTll998l : ISmith et aLlbOOTl ). For Gaussian initial 
conditions, this has been widely tested against numerica l 
simulations (e.g. ISheth et al.1 l200ll : ISeliak fc WarrenI |2004| : 
[Tinker et aLll2005l v 

In general, the halo bias can be quantified using either 
the power spectrum of the halo density field, Phh, or the 
cross-spectrum between the halo and the underlying matter 
density field, Phm- In the two cases the bias reads 



6hh(fe,M,z) 



6hm(fc,M,z) 



Phh(fc,M,z) 
PiKz) ' 

Ph^(fc,M,z) 
P{k,z) ' 



(12) 



(13) 
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Figure 7. Ratio between tlie matter power spectra of non- 
Gaussian and Gaussian simulations at redshift z = 1 (bottom) 
and 2 = (top). Data arc extracted from the N-body simulations 
in our main series where identical random phases have been used 
to generate (j> for all values of /nl- 



where P{k, z) is the matter power spectrum. If the bias due 
to halo formation is local and deterministic then 6hh = &hm 
apart from measurement errors. However, in the presence of 
a stochastic component that does not correlate with the den- 
sity field fehh ^ &hm- In practice, however, the measurement 
of all power spectra is affected to some level by shot noise 
due to the discrete nature of dark-matter haloes and N-body 
particles. If the distribution of the tracers can be approxi- 
mated as the Poisson sampling of an ideal density field, then 
the measured power spectrum corresponds to tha t of the 
unde rlying field plus the mean volume per particle l|Peeblesl 
Il98d ). Discreteness effects are thus expected to be negligible 
for P and Phm due to the large number density of particles 
in the simulations. On the other hand, massive haloes are 
rare and, being extended objects, cannot be m odelled as the 
Poisson sampling of a continuous dist ribution (|Mo &: White! 
1 19961 . iMaghocchetti fc Porciaiiil |2003| . Porciani in prepara- 
tion). It is not clear then how to correct for the d iscreteness 
effect in their power spectrum (|Smith et al.1120071 ) . For these 
reasons we use &hm in our analysis and we adopt tihh (with- 
out performing any discreteness correction) only to verify 
the results (see Figure IS)). 



5.1 Halo bias from Gaussian initial conditions 

It is well known that the halo bias factor from Gaussian ini- 
tial conditions is approximately scale-independent for small 
values of the wavenumber fc. We will refer to this asymptotic 
value on large scales as the "linear bias" and denote it by 
bo- Similarly to the halo mass function, when expressed in 
terms of o"~^, the linear bias assumes a universal form which. 
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Figure 8. The halo bias from the halo-halo power spectrum (with 
no discreteness corrections) is plotted against the halo bias from 
the halo-matter cross spectrum. Whenever the density of haloes 
is high enough, the two estimates are very close showing that 
little stochasticity between mass and halo overdensities is present 
on the scales of interest (indicated in h Mpc~^ in the label). 
The excess in fejih for rare, massive haloes is likely due to shot 
noise. Note that large positive values of /nl correspond to more 
massive haloes and thus allow more accurate measures of high 
bias parameters. 



within a given accuracy, is independ ent of redshift and just 
weakly dep endent on cosmology (e.g. lSheth fc Tormen|[l999l : 
ISeliak fc Warrcn..2004i ). 

We measure the linear bias for the haloes in our simu- 
lations as follows. We first determine the functions bhh and 
bhui by directly applying equations (|12|l and (|13|l . Within the 
statistical uncertainties, both functions approach asymptot- 
ically to a constant on large scales [k < 0.05 /i Mpc~^). We 
use the average of the bias function measured in the range 
0.01 < k < 0.05 h Mpc~^ (4 fc-bins) as our estimate of the 
linear bias. The standard error of the mean is used to quan- 
tify the corresponding statistical uncertainty0 

In Figure [9] we show the linear bias obtained from Run 
1.0 (triangles). Run 2.0 (squares) and Run 3.0 (circles) as 
a function of . Simulation data from snapshots between 
z — and 2 = 2 are compared with the commonly used 
parameterisations listed in Table [7] Our re sults are in good 
agreement with th e fit bvlSheth et al.l (l200ll ) for large masses 
and with that by [Tinker et al. ( 2003) for smaller masses. 
Note that by combining simulation boxes we are able to 
explore a larger interval of a^^ than previous studies. 

Given that no existing model for the linear bias ac- 
curately reproduces our results over the entire mass range 
spanned by the simulations, we decided to derive a new fit- 
ting formula. In particular, we parameterised the outcome 



^ Consistently, in what follows, we use the rms value as the error 
on b(fc, M, 2, /nl = 0). For /nl Oi we assume that the relative 
error is the same as in the Gaussian case. 
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Table 7. Commonly used parameterisations for the linear bias arising from Gaussian initial conditions. 



Acronym 



Reference 



Functional form 



Parameters and Variables 



MW 



ST 



sw 

T 



Mo fc White (19961 
Sheth et al. ('20011 



^MW = 1 + ■ 



fegrp = 1 + 



+ y/ab I a 



+ 6{l-c){l-c/2) 



Seliak fc Warren (20041 bsw = 0.53 + 0.39{x°'^^) + j^^^f^ 



Tinker et al. (20051 



+ 5 X lO-'^x^ 



Sc = 1.686 
Sc = 1.686 

a = 0.707,6 = 0.5, c = 0.6 

x= jf- with a{M,) = 1.686 
1.686, a = 0.707, b = 0.35, c = 0.80 
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Figure 9. Linear halo bias from the Gaussian simulations Run 
1.0 (triangles). Run 2.0 (squares), and Run 3.0 (circles), as a 
function of a~^. Only bins containing more than 1000 haloes are 
shown. The solid lines correspond to the functions listed in Table 
[T] as indicated by the labels. The four hex agons correspond to 
the data at z = 10 bv lGohn fc White! bOOSll . The vertical dotted 
lines indicate the maximum and minimum con sidered by 

[Tinker et al] | |2005|) (red) and ISeliak fc WarrenI l[20oi) (cyan, in 
this case the minimum cr"^ coincides with the frame of the figure). 
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Figure 10. Scale-dependent halo bias arising from non-Gaussian 
initial conditions. Results are shown in terms of the ratio between 
the bias functions measured from a simulation with a given /jml 
and with /nl = at fixed halo mass (indicated by the label in 
units of /i~^Mq). Note that in the Gaussian case the bias keeps 
nearly constant for k < 0.05 h Mpc~^. 



and used minimisation to find 



Bo 

B2 



0.647 ±0.010 
-0.540 ± 0.028 
1.614 ±0.019 



(15) 



• i„+-„rio Qo This fit (which reproduces the numerical data with great 

accuracy in the range -1.1 < Incr"' < 0.8) should be con- 
sidered as the linear bias naturally associated with the mass 
bo = Bo + Bi + B2 <y~^ , (14) function given in equations (jl} and 
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5.2 Halo bias from non-Gaussian initial conditions 

Recent analytical models, have suggested that the halo bias 
arising from non-Gaussian initial conditions of the local type 
does not tend to a constant on large scales. Rather, the 
deviation from the Gaussian case should follow 



Ah = 6(fc,M,2,/NL)-fo(fc,Af,2,/NL=0) = 

5c g{oo) Hi n„ 



= 3/nl[&o(M,2)~1]. 



, (16) 



D{z) g(0) c2 k^T{k) 

where 6c = 1.686, c/Hq = 2997.9/1"^ Mpc is the Hubble 
radius, T{k) is the matter transfer function, and D{z) is 
the linear growth fac tor of matter perturbations normalised 



2008 




2008 


) 



Slosar et al]|2008l : lAfshordi fc Tollevll2008l: iMcDonald 



have indeed shown that the halo bias is scale dependent 
even for small values of k in non- Gaussian cosmologies (with 
|/nl| = 100,500) and found qualitative agreement with 
equation (|16p . In Figure [Till we show how the bias depends 
on scale in our simulations which also consider smaller val- 
ues of |/nl|. Our results confirm the presence of a strongly 
scale-dependent bias. Larger values of |/nl| correspond to 
a more marked scale dependence. Note, however, that for 
k > 0.05/1 Mpc"^ the non-Gaussian deviation Ab changes 
sign. On these scales, the halo-matter and halo-halo spectra 
emerging from non-Gaussian perturbations has actually less 
power than in the Gaussian case. The opposite happens with 
the matter power spectrum (even to a larger degree) and the 
net effect is a negative Ab. This result implies that equation 
(|16p can only hold asymptotically on very large scales. This 
is not suprising if interpreted within the peak-background- 
split formalism where the large-scale bias is linked to the 
first derivative of the mass function with respect to . In 
the non-Gaussian case the bias is composed of two parts, a 
scale- independent term and the correction given in equation 
(|16p. Since the halo mass function changes shape when /nl 
is varied, also the constant bias should depend on /nl for 
a fixed halo mass. Increasing /nl corresponds to a larger 
abundance of massive haloes and to a slightly smaller con- 
stant bias with respect to the Gaussian case. Likely, this is 
what makes the Afe in the simulations negative for positive 
/nl. To proceed with a detailed analysis of our simulations, 
we find it convenient to rewrite equation (|16|l as 



Ab = /nl (6o - 1) 



a{k, z) 



(17) 



where F — 2i ScVlui Hq / c? and a(fc,2) = 
k^T{k)D{z)g{Q)/g{oo). In Figure [11] we test the scal- 
ing of Ah with redshift, linear bias and wavenumber for 
/nl = +750 (where we have the best signal-to-noise ratio 
at high halo mass). Similar results are obtained with 
different values of /nl. The quantity shown is A6a/F which 
should correspond to /nl (6o ^ 1) if the analytical model 
provides a good description of the data. This quantity is 
indicated by a dashed line. The following two trends clearly 
emerge from the data. For small values of k, the model 



^ The facto r g(oo )/g(0) is needed since [Dalai et al. ] l l2008t) and 
[Slosar et al.l bOOSi ) normalise the growth factor D(2) to be (1 -|- 
z)~^ during matter domination. 



overestimates the data by 20-70 per cent increasing with 
&o and independently of z. On smaller scales, discrepancies 
become more and more severe. At ~ 0.05 h Mpc~^, 
the model is systematically a factor of 5 higher than the 
data. The fc-dependence of A6 is therefore different than in 
equation (|16|) . 



The data also drop a hint that, for k > 0.01 h Mpc~^, 
the scaling with &o — 1 might only persist up to a maximum 
value of 6o, 6o,max. For bo > bo, max it appears that the 
values of Afe are always smaller than expected from the 
extrapolation of the trend bo — 1 determined at smaller 
fee. The value of bo, max seems to depend both on redshift 
and wavenumber and roughly corresponds to constant 
halo mass for a given k. However, uncertainties in Ab at 
these high masses become very large and it is difficult to 
judge how robust the existence of bo, max really is. We note 
anyway that when we tried to fit data at different redshifts 
(for a given /nl and k > 0.015 h Mpc~^) by adding a 
variable normalisation constant in front of equation (|16p . 
we systematically obtained significantly different fits (at a 
confidence level of 2.5 a) at different redshifts. This trend 
disappears when only the lowest values of bo are considered 
at each redshift for the fit. 

Data from simulations with all the considered values of 
/nl are shown with different symbols and colors in Figure 
1121 Each panel refers to a particular wavenumber bin (in- 
dicated by the label in units of h Mpc^^). The model in 
equation (|16|) is again indicated by a dashed line. Note that, 
in most cases, it substantially deviates from the simulation 
data. In particular, Ab measured from the simulations shows 
a much stronger fc-dependence than the analytical formula, 
as already seen in Figure 1111 In general, the overall ampli- 
tude of Ab drops by an extra factor of ~ 3 with respect to 
fc^ T{k) when moving from k ~ 0.01 h Mpc"^ to fc 0.05 h 
Mpc~^ independently of bo and /nl- Also, Ab does not seem 
to scale linearly with /nl while its linear dependence on 
bo — 1 appears to be solid, at least for bo < bp, max. We thus 
introduce a correcting factor /3(/NL,fc) defined by 

/3(/NL,fc)/NL(bo-l)^^^ (18) 



Ab: 



a{k, z) 



and we measure it by fitting the simulation data for 
b(fc,M, 2, /nl) and b(fc,Af, 2,0) at constant values of /nl 
and k. We use an effective variance weighted least squares 
method to simultaneously account for errorbars on both bias 
parameters. The best-fitting values are reported in Table [8] 
and can be used to compute the function /3 by interpolation. 
The final expression for Ab, corrected with the /3 factor, is 
shown in Figure [12] with solid lines. 



Data in Table[8]have an amazing regularity. Apart from 
a normalisation constant, each column (row) shows the same 
linear trend with k (/nl). This suggests that, within the 
explored parameter range (0.01 < k < 0.05 Mpc~^ and 
-80 ^ /nl =S 750), 

P{k, f^i^) = l3o (1 - /3i /nl) (1 -P2k) . (19) 

We thus use this equation to fit the original data for the halo 
bias from Gaussian and non-Gaussian initial conditions and 
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Figure 11. The non-Gaussian bias correction Ab as a function of the linear Gaussian bias feo for /nl = 750. Data from the simulations 
are indicated by symbols with errorbars and correspond to different redshifts as indicated in the label. The dashed line marks the 
prediction of the model in equation HIGH . The wavenumber in the labels is given in units of h Mpc~^. Horizontal errorbars are not drawn 
to improve readability. 



Table 8. Best-fitting value and Icr uncertainties for the multiplicative correction /9(A:, /nl). The first set of data corresponds to the 
fc-intcrval where the Gaussian bias is constant. 



fc(/iMpc-i) 


/3(fe,-80) 


/3(fc, -27) 


/3(fc, +27) 


/3(fc,+80) 


/3(fc,+250) 


/3(fc,+500) 


/3(fc,+750) 


0.0117 


0.97 ± 0.04 


0.88 ±0.11 


0.83 ±0.13 


0.83 ±0.05 


0.71 ± 0.02 


0.66 ±0.01 


0.60 ± 0.01 


0.0191 


0.77 ± 0.07 


0.78 ±0.21 


0.69 ±0.25 


0.70 ±0.10 


0.63 ±0.03 


0.56 ±0.01 


0.51 ± 0.01 


0.0303 


0.65 ±0.14 


0.58 ±0.27 


0.59 ±0.27 


0.54 ±0.16 


0.50 ±0.06 


0.42 ± 0.03 


0.37 ± 0.02 


0.0494 


0.35 ±0.32 


0.45 ±0.19 


0.40 ±0.26 


0.31 ± 0.33 


0.29 ±0.09 


0.25 ±0.05 


0.20 ± 0.02 


0.0804 


0.01 ± 0.23 


0.20 ±0.21 


-0.10 ±0.43 


-0.06 ±0.22 


-0.08 ±0.18 


-0.10 ± 0.09 


-0.11 ± 0.04 
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find 

f3o = 1.029 ± 0.027 , 

/3i = (4.25 ± 0.33) X 10"* , (20) 
/32 = 14.8 ± 0.5 h'^ Mpc , 

at the 68.3 per cent confidence level. Note that we computed 
the power spectra in finite-sized bins of the wavenumber, so 
that there is some degree of ambiguity in associating the re- 
sults with a given value of k. Unfortunately the choice plays 
a role in determining /3 as a is a steep function of k on the 
scales of interest. In Tableland in equation (|20|) . we have 
used the arithmetic mean of the wavenumbers contributing 
to a given bin. If one instead uses the logarithmic center 
of the bin, /3o is slightly reduced with a best-fitting value 
of 0.897 ± 0.024. The parameters /3i and 02 are unaltered. 
Therefore, a systematic contribution ~ 0.1 should be added 
to the error budget of Pq. 



equations (|16|l and 1)18^ assume that the Gaussian bias 
bo is constant with k but this is only approximately true in 
the simulations. The fit in equation (|20|l . the Table|S]and the 
Figures (lllf) and (|12|) have been obtained by identifying bo 
with the actual bias measured in the Gaussian simulation at 
each wavenumber. If, instead, the estimate for bo introduced 
in Section 5.1 is used, one gets f3o — 0.970 ± 0.027, /3i — 
(4.13 ± 0.33) X lO""* and 02 = 13.8 ± 0.7 h'^ Mpc, slightly 
improving the goodness of the fit. 

The quadratic dependence of Ab on /nl is rather sur- 
prising as it cannot be straightforwardly derived from the 
simple models listed above. It might possibly arise from 
higher-order terms which have been neglected in the expan- 
sion that leads to equation (|16fl . Anyway, it is clearly present 
in the simulations as it can be seen by looking at the varia- 
tion of /3 along a given row in Table |5] Within the range of 
/nl of physical interest, the effect is rather small; the coeffi- 
cient Pi only corresponds to a few percent correction. Note 
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that a quadratic term breaks the symmetry in the ampli- 
tude of A6 between non-hnearity parameters with opposite 
sign and identical absolute value. It is hard to directly test 
this against our simulations as we just have two runs with 
/nl < and both of them correspond to rather small |/nl| 
where the uncertainties in /3 are large. An alternative ex- 
planation for a non- vanishing /3i could be that it artificially 
derives from imposing a linear relation in bo — 1 to data that 
do not scale linearly for 6o > feo.max. Indeed, just using dat- 
apoints with small values of 6o we derive bigger values of /3 
for large /nl (more or less in line with /3i =0). Therefore, 
what is robust is that at least one of the scalings with bo 
or with /nl is incorrect in equation (|16|) . We found that a 
scaling proportional to 70 (1 -I- 71 log 60) (with 70 and 71 two 
adjustable parameters) does slightly better (in terms of re- 
duced x^) than Po (60 - 1), at least for k > 0.014 ft Mpc"^ 
Howe ver, since the scaling with 60 — 1 has a sou nd theoretical 
basis (|Mo fc Whitelll99g : ICatelan et al.lll998l ) we preferred 
to quote our results as in equation (|19|l . From the statisti- 
cal point of view, the parameters pO|) provide an accept- 
able description of the simulation data to high confidence 
for all values of 60. However, they are particularly accurate 
for 60 > 2 — 2.5, while /3i ~ (with the same /3o and (52) has 
to preferred for smaller values of 60. 

The linear correction in k should be thought of as the 
first-order term of a series expansion in the wavenumber. 
We attempted to determine the corresponding quadratic 
term by considering larger values of k in the fit (one bin 
more, up to fc = 0.0962 ft Mpc~^). However, values of A& 
become small compared with the numerical errors and we 
found that the quadratic parameter is badly constrained 
by the data (/Sa = 34 ± 34 ft"^ Mpc^) while the other 
parameters remain nearly unchanged (and get larger uncer- 
tainties). Also note that the Gaussian bias starts to depart 
from 60 at fc > 0.05 ft Mpc^^ and it is not clear whether 
equation (|16p should still be expected to hold in this regime. 

iDalal et all (|2008l ') derived an expression for Afe which 
coincides with equation (|16|l but does not include the linear 
transfer function. Theoretically, this is hard to understand, 
as non-Gaussianity is generated well before matter-radiation 
equality and one should account for the linear evolution 
of density perturbations. Anyway, due to the different fc- 
dependence, their expression for A6 provides a better fit to 
the simulation data than equation (|16|l when both models 
are allowed to vary in amplitude with a tunable free param- 
eter0 None of them, however, provides such an accurate fit 
to the data as our equations p9|) and (|20[) . which improve 
the by at least a factor of 1.7. 



6 DISCUSSION 

Slosar et al. (2008) have used equation (|16|) to constrain 
/nl by considering measures of the clustering amplitude of 
luminous red galaxies (LRGs) and quasars from the Sloan 
Digital Sky Survey. Combining all datasets, they found 
—29 < /nl < +70 at 95 per cent confidence. How would 

* The b est-fitting value f or this coefficient reads 0.92 for the 
model of lDalal et all ||2008D and 0.58 for the model of lSlosar et al.l 
1I2OO8I ') 



this result change based on our simulations? Disentangling 
the different contributions, the strongest constraints to /nl 
in Slosar et al. (2008) come from the angular power spec- 
trum of quasars with photometric redshifts in the range 
1.45 < z < 2.00 and a mean bias of ~ 2.7. Weaker lim- 
its are also contributed from the power spectrum of spec- 
troscopic LRGs and the angular spectrum of photometric 
LRGs. (with a bias of ~ 2 at 2 ~ 0.5). Figure HH) and Ta- 
ble[S]suggest that at the scales of interest (0.01 < fc < 0.05 ft 
Mpc"^) the model given in equation (|16|) tends to over- 
estimate the scale-dependent bias seen in the simulations. 
Therefore, to match an observed A6, a larger value of |/nl| 
is required than predicted by the analytic model. When ap- 
plied to the data by Slosar et al. (2008), our correction 
should thus give somewhat looser limits on /nl. Because 
of the strong fc-dependence of the function (3 it is impossible 
to give more precise estimates without fitting the power- 
spectrum data. Note, however, that a steeper fc-dependence 
potentially makes determinations of /nl even more compet- 
itive with respect to studies of CMB anisotropies. 



7 SUMMARY 

We use a series of high-resolution N-body simulations 
to study the mass function and the clustering proper- 
ties of dark-matter haloes arising from Gaussian and 
non-Gaussian initial conditions. In particular, we follow 
the formation of structure in a universe characterised 
by the best-fitting parameters from the third- and fifth- 
year WMAP data releases. We consider non-Gaussianity 
of the local type and we use eight different values of /nl 
(-80, -27, 0, +27, +80, +250, +500, +750) enclosing the pa- 
rameter space currently allowed by studies of the cosmic mi- 
crowave background. Our main results can be summarised 
as follows. 

(i) The mass function of dark-matter haloes varies with 
/nl. Larger values of the non-linearity parameter corre- 
spond to higher abundances of the most massive haloes. 
Analytical models ba s ed on the Press-Sc hechter method 
iMatarrese et all l200d : ILo Verde et alll2008h are compatible 
with our simulated results for the ratio of the Gaussian and 
non-Gaussian mass functions only if the critical threshold 
for halo collapse is lowered to 5c ~ 1.5. An accurate fit of 
the Gaussian dn/dM is necessary to derive the non-Gaussian 
mass function from the aforementioned ratio. 

(ii) We find that, in per fect analogy with the Gaus- 
sian case (| Jenkins et al.|[2001^ . the halo mass function as- 
sumes an approximately universal form. This means that, 
when expressed in terms of suitable variables, its depen- 
dence on redshift and cosmology is erased to good precision 
(nearly 10 per cent). We parameterise the /NL-dependence 
of the universal mass function and provide an accurate fit 
for its high-mass end. For —80 ^ /nl ^ 250 and for masses 
M > lO^'^ft~^M0, the best-fitting parameters for the non- 
Gaussian halo mass function in equation Q are given in 
equation (O and Tabled This fit reproduces the mass func- 
tion of friend-of-friends haloes with an accuracy of 5 per 
cent on top of a systematic contribution (up to 10 per cent) 
due to the non perfect universality. For applications requir- 
ing higher precision, an additional formula is provided: for 
-80 ^ /nl s; 80 and ^ 2 ^ 0.5 the fit in equations (fTO)) . 
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((U, and (lllf) has to be preferred to the universal fit. On the 
other hand, for higher values of |/nl| and for higher red- 
shifts, the universal fit gives a better and more economic (in 
terms of parameters) description of the data. In the Gaus- 
sian case, we extend the fit to a larger interval of halo masses 
[M > 2.4 X 10^° /I'^Mo) by combining simulations with dif- 
ferent box sizes: - see equations Q and ([5]). Our fitting 
function provides a precious tool to forecast constraints on 
/nl from future surveys and to analyze current datasets. 

(iii) The matter power-spectrum in non-Gaussian cos- 
mologies departs from the Gaussian one already on very 
large scales. For values of /nl within the current CMB con- 
straints these scale-dependent deviations can be as high 
as two per cent at fc = 0.3 h Mpc~^ and increase with 
wavenumber. The discrepancy is systematic: models with 
positive /nl have more large-scale power than the Gaussian 
case and models with negative /nl have less. This warns 
against the widespread habit of using the Gaussian matter 
power spectrum to determine non- Gaussian bias parameters 
when high-precision is required. It also suggests that pri- 
mordial non-Gaussianity modifies the shape and amplitude 
of the baryonic-oscillation feature in the two-point statistics 
and the convergence power spectrum in weak-lensing stud- 
ies. 

(iv) We present an accurate fitting formula for the lin- 
ear bias of dark matter haloes arising from Gaussian initial 
conditions extending previous work to larger mass intervals. 
This, together with the mass function fit mentioned above, 
can be used to constrain parameters of halo-occupation 
models from clustering data. 

(v) Finally, using the halo-matter cross spectrum, we 
confirm the strong fc-dependence of the halo bias on large 
scales (fc < 0.05 h Mpc~^) which was already detected by 
iDalal et al. I (|2008t ). However, we show that commonly used 
parameterisations based on the peak-background split over- 
estimate the effect for k > 0.01 ft Mpc^^. The discrepancy 
increases with the wavenumber and at fc > 0.05 /i Mpc~^ 
A6 in the simulations changes sign with respect to the 
models. On top of this, the analytic model for the scale- 
dependent part of the bias requires corrections which de- 
pend on the non-linearity parameter, the wavenumber and, 
possibly, also on redshift and clustering strength, equations 
H18p and p9p with the best-fitting parameters listed in (|20p 
provide a fitting formula which accurately reproduces the 
outcome of the simulations for 0.01 < k < 0.05 h Mpc~^ and 
^80 ^ /nl ^ 750. This fit should be employed to constrain 
/nl from future clustering data at low and high redshift. 
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APPENDIX A: INITIAL CONDITIONS AND 
ZEL'DOVICH TRANSIENTS 

The initial positions and velocities of the particles in our N- 
body simulations have been generated using the Zel'dovich 
approximation. This method introduces long-lasting artifi- 
cial transients in the growth of perturbations which might 
alter the halo m ass function at the high-mass end even at 
very late epochs (jCrocce et al.|[2006l l. It is therefore impor- 
tant to start the simulation at a sufficiently high redshift 
to ensure that all transients have decayed within the cos- 
mic time at which the simulation output is used for science 
applications. Alternatively, less stringent requirements on 
the initial redshift are necessary if one uses second order 
Lagrangian perturbation theory to displace particles at the 
initial time (jCrocce et al.|[2006f ). 

The simpler Zel'dovich approximation (which only re- 
quires the calculation of the gravitational potential) is much 
more widespread. In this case, a few authors have inves- 
tiga ted how to compu te the optimal starting redshift (see 
e.g. iLukic et al.] |2007^ as well as quantified the effects of 
the i nitial redshift on the halo mass function (| Tinker et al.l 
l2008h and on the i nternal properties of dark matter haloes 
jKnebe et alll2009l '). These studies show that simulations of 
the concordance cosmology (and gaussian initial conditions) 
with initial redshifts of 35 < Zstart < 60 (depending on the 
simulations specifications) have converged to the correct so- 
lution by 2 = (at least for halo masses M < 10^'* h^^ AIq). 
Even though our Zstart is in the right ball park, it is impor- 
tant to test that our results are robust against changing 
initial redshift. We thus decided to perform the following 
simple test: we re-simulated Runl.O and Runl.750 of our 
main Series using aatart=99 instead of 2atart=50 and com- 
pared the halo mass functions of the two simulations as a 
function of redshift (see Figure Al). In good agreement with 
[Tinker et al.l (|2008l l . for /nl = we find that discrepancies 
are smaller than 10 per cent at z ~ 2 and that the (pos- 
sible) effects of Zel'dovich transients are completely erased 
by z ~ 1. Our results show that this holds true also for 
relatively large values of /nl (see the right panel in Figure 
Al where /nl = -1-750): even though non-linearities in the 
initial conditions are slightly enhanced with respect to the 
Gaussian case, at any given redshift the corresponding den- 
sity field is also more evolved and the effect of the Zel'dovich 
transients on the halo mass function are again erased by 
z ~ 2. For this reason this paper only uses data from snap- 
shots at 2 ^5 1.6 which are accurate to better than 5 per 
cent. 
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Figure Al. Effect of the initial redsliift ^start on the halo mass- function, for /nl = (left panel) and /nl = -1-750 (right panel). 



